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l/") ' Motivated by recent experiments we have carried out an Event Driven computer simulation of 

a diluted binary mixture of granular particles vertically vibrated in the presence of gravity. The 
simulations not only confirm that the kinetic energies of the two species are not equally distributed, 
as predicted by various theoretical models, but also seem to reproduce rather well the density and 
^ , temperature profiles measured experimentally. Rotational degrees of freedom do not seem to play 

. any important qualitative role. Instead, simulation shows the onset of a clustering instability along 

the horizontal direction. At the interior of the cluster we observe a secondary instability with respect 
to the perfect mixing situation, so that segregation of species is observed within the cluster. 
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I. INTRODUCTION 

a 

' The present keen interest in the dynamical properties of granular materials is motivated both by the challenge of 
understanding the complex processes involved and by the important practical applications in engineering, industry 
and technology fl . These materials are peculiar in many respects and display several intriguing phenomena such as 
clustering , shear instability j^, Q and lack of energy equipartition, which make their behavior different from ordinary 
molecular fluids. The dissipation of kinetic energy during the inelastic collisions makes them special. The main 
motivation of the present paper stems from two recent experiments @] , p] which demonstrated that when a mixture 
constituted by two different species of grains is vibrated, each component attains its own "granular temperature", 
i.e. the average kinetic energy per particle does not take on the universal value fKT, where / is the number of 
degrees of freedom and K a constant, as it occurs in molecular gases. On the contrary, one observes that the ratio 
Z1/X2 varies with concentration, inelasticity parameters, particle sizes, masses and driving mechanism. Even in 
the absence of energy injection, the inelastic gas cools, but one observes that the temperature ratio asymptotically 
^ remains constant. On the other hand, it is understood that while the only relevant hydrodynamic field is the global 
C^} temperature T = xT\ + (1 — x)T2, transport properties depend on that ratio 0. 

"jS ■ These are of course manifestations of the non-equilibrium nature of Granular systems, which can only be maintained 
stationary by a continuous energy feeding to compensate the energy losses due to the inelastic collisions and to friction. 

A theoretical understanding of such a behavior of granular mixtures has been achieved in the case of homogeneous 
driving mechanisms by means of a combination of models and approximations including the pseudo-Maxwell inelastic 
gas and the Inelastic Hard Sphere model treated by means of the Boltzmann-Enskog equation. Both models have 
been studied analytically and numerically in the free cooling || ^| and in the driven case (l0| [n], [l2| Apart 
from the studies of ref JlJ, |l3| none of these investigations considered the role played by the gravitational field, by 
the strongly inhomogeneous boundary conditions employed in the experiments of refs. ^, |(|, by the roughness of the 
grains and by their rotational degrees of freedom. Moreover, there is still an open debate about the "best" energy 
feeding mechanism. Whereas theoreticians seem to favour a uniform thermal gaussian bath, because it lends itself 
to a great deal of analytical work, a numerical computer experiment can test directly driving mechanisms which are 
closer to those employed in a laboratory. The structure of the paper is the following: in section 2 we illustrate briefly 
the model, leaving the technical details to the appendix; in section 3 we discuss the results for the geometry of ref. [|| . 
In sec. 4 we consider a different aspect ratio, namely a longer box, where gravity plays a more relevant role. Finally 
in sec. 5 we present our conclusions. 
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II. MODEL SYSTEM 

We decided to remain as close as possible to a commonly employed experimental set-up, by constraining the grains 
to move on a vertical rectangular domain of dimensions L x x L y . The gravitational force acts along the negative y 
direction. The grains are assumed to be spherical and free to rotate about an axis normal to the xy plane. They 
receive energy by colliding with the horizontal walls, harmonically vibrating at frequency v. The side walls instead 
are immobile and were chosen either smooth or rough according to the numerical experiment. When side walls are 
considered to be rough, they are assumed to have the same friction coefficient /i as the particles. 

The collisional model adopted in the present paper corresponds to the one proposed by Walton Q. ft conserves 
both the linear and the angular momentum of a colliding pair, but allows energy to be dissipated by means of a 
normal restitution coefficient and a friction coefficient fi. The collision rule (given in detail in the appendix) takes 
into account a reduction of normal relative velocity of the two particles (V n ), a reduction of total tangential relative 
velocity (V r ) and an exchange of energy between those two degrees of freedom. The reduction of normal relative 
velocity is modeled by means of a non constant restitution coefficients CKy S [0, 1], whose dependence by the relative 
velocity is of the form: 




for V n < vq 



{ Tij for V n > v 

where i and j are the numbers indicating the species of the colliding particles, ry are constants related to the three 
types of colliding pairs, t>o ~ \> r gd, where d is the average diameter of the particles and g the gravitational acceleration 

Simulated collisions are of two types: with sliding or sticking point of contact. When the following condition is 
satisfied (high relative tangential velocity), the collision happens in a sliding fashion, otherwise it is sticking: 

^>^p/*(l + «tf) (1) 

where I is the dimensionless moment of inertia (equal to 2/5 for spheres), while fi is a static friction coefficient 
characterizing the surface roughness of particles, assumed equal to the dynamical friction coefficient. The full dynamics 
consists of interparticle collisions, and wall-particle collisions. The trajectories between collision events are parabolic 
arcs due to the presence of the gravitational field. 

An efficient Event Driven (ED) simulation code was employed to evolve the system |f(J]. 

The two species were chosen to be spheres of equal diameters d = 0.16cm, and unequal masses mi = 1 .58 ■ f 0~ 2 g and 
m-2 = 5.21 • 10~ 3 <?, respectively. The driving frequency was set to 50 Hz, the vibration amplitude A — 3.5 diameters 
so that the corresponding dimensionless acceleration T = Au 2 / g — 56. 

All the averages quantities reported in the following have been obtained by employing equally spaced data points 
separated by time intervals At = fCP 1 s in order to assure statistical independence of measures. We performed a 
number n = 1.5 • 10 5 of vibration cycles. 

Barrat and Trizac Jl3|] have recently considered one of the systems studied in the present paper. However, our 
treatment presents some differences: 

• The bottom and top walls of |p"3| move in a sawtooth manner with a negligible excursion so that their positions 
are considered fixed, while our walls move sinusoidally with a non-negligible amplitude. 

• The walls are not smooth in our treatment, but have a friction coefficient /i > 0. 

• The collisional models of our treatment and that of Barrat-Trizac are different (they do not distinguish between 
sticking and sliding collisions). 

• We take into account gravity. 



III. RESULTS FOR A SHORT SYSTEM 

The box dimensions were L x = 48 d and L y = 32 d. Simulation runs were carried out using Ni = 150 grains of each 
species. The static friction coefficient has been always chosen as /i = 0.1. The stationary state is determined by the 
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FIG. 1: Typical snapshot of the system described in III. Open circles indicate particles of species 1, black circles particles of 
species 2. 




FIG. 2: Top panel: Translational temperature profiles for species 1 (circles) and species 2 (squares) for a mass ratio mz/mi = 
3.03, expressed in fiJ (left scale), and temperature ratio (diamonds) T2/T1 (right scale). The vertical position is measured in 
particle diameters (d) relative to the geometric center of the cell. Bottom panel: rotational temperature profiles for species 1 
(circles) and species 2 (squares) and temperature ratio T2/T1 (diamonds). 



balance between the energy input provided by the vibrating walls at a frequency of 50 Hz and the dissipation due to 
inelastic collisions. The typical collision frequencies are of the order of v\ ~ 580 Hz and V2 ~ 850 Hz for the heavy 
and light balls, respectively. A typical microscopic configuration of the system is shown in fig.|l|. 
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FIG. 3: Same as in fig. 1, but with different restitution coefficients. 

A. Temperature profiles 

The restitution coefficients were first set to r\\ — 0.93 for 1-1 collisions, r^2 — 0.7 for 2-2, and r%2 — 0.8 for 
1-2 collisions. This means that the more massive particles are also the more elastic ones. In figure || we show the 
partial translational temperature profiles for the two species, and observe that close to the vertical boundaries the two 
temperatures are essentially determined by the energy injected by the vibrating walls. Indeed, interparticle collisions 
are rare within this region, and play no signifcant role because of the low local density (see figure |J). In addition, 
grains 1 and 2 impinging with the same speed on the mobile wall bounce with the same velocity (V oc Ad), hence the 
local value of the temperature ratio, 7 = T2/X1 near the vibrating walls, turns out to be approximately 7 ~ mi/m^, 
as shown in figure |[ 

On the other hand, the temperature drops as the distance from the walls increases, while the ratio 7 grows up to 
a plateau value, indicating that collisions tend to cool the mixture and render the two partial temperatures closer. 
Figure ^| clearly displays the breakdown of the kinetic energy equipartition already noticed in previous experimental 
and theoretical studies. 

We also measured the rotational temperature profiles, shown in figure |^. We observe for rotational temperatures the 
same kind of equipartition breakdown that holds for translational degrees of freedom. Moreover, the absolute values of 
rotational and translational temperatures are quite different, as already reported by Luding flSfl for a one-component, 
vibrated granulate. On the other hand, the ratio of the two rotational temperature profiles seems to be quite close to 
that of the translational temperature profiles. 

In figure^, we changed the restitution coefficients and set r\\ = 0.7, T22 — 0.93 for 2-2, and ri2 = 0.81. Therefore 
now the more massive particles are the more inelastic. We observe that, whereas the temperature profiles near the 
vibrating planes are nearly unchanged, because collisions are rare, the value of the temperature of the heavier species 
is lower and the temperature ratio is closer to 1. In this case the larger inelasticity of the heavier particles competes 
with the mass asymmetry which instead tends to make T2/T1 smaller. 

B. Area fraction profiles 

Due to the large value of the parameter T, the partial density profiles tend to be rather symmetric, with a maximum 
near the center. We also notice small differences in the density profiles, which reveal that the heavier species has 
higher concentration close to the center of the cell, while the lighter species is more spread. 

Comparing the present results for the area fraction profiles with those recently obtained by employing a Direct 
Simulation Monte Carlo technique JL2| , we notice that the agreement is only qualitative, whereas the temperature 
ratios are in significantly better agreement. 

We verified that removing tangential friction (coupling the translational and rotational degrees of freedom) does 
not change significantly the above scenario. In fact, the translational temperature and density profiles of the cases 
fx = and /1 = 0.1 show only small quantitative differences, more pronounced close to the horizontal walls. 
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FIG. 4: Area fraction profiles m and ri2. Control parameters are the same as in fig. 1. 
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FIG. 6: Same as in figure El but with different restitution coefficients. 
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FIG. 7: Rescaled velocity distribution functions < c > P(c) = ,f(c/ < c >) for translational velocities and angular velocities 
of both species. Set A corresponds to parameters of fig. 2, set B to parameters of fig. 3. The independent variables (v, u) 
were rescaled by their mean squared value (no, uio). On rescaling, the distributions collapse nicely onto each other. The dashed 
curve is the fitting law discussed in the text, the continuous line is the Gaussian plotted as a reference for the eye. 



C. Transversal profiles 



We also studied the density and temperature profiles along the horizontal direction. To the best of out knowledge, 
no such measure has been reported in experimental works. In figures |5| and ^ we observe that temperature profiles 
vanish close to side walls, while density profiles display their maxima in the same region. In order to gain further 
insight, we analyzed a sequence of snaphsots of the dynamics, and observed that the system bears a denser cloud of 
grains in the vicinity of one of the side walls. Such a configuration was maintained over an interval of time much 
longer than the vibration period 2tt/lu. 

The cluster eventually "evaporates" to form again close to a randomly selected side wall. Over several periods 
of oscillation of the cell, we noticed an effective horizontal simmetry breaking, i.e. the number of particles in the 
right-hand and left-hand sides of the cell were rather different. The system is unstable with respect to horizontal 
density fluctuations and clusterizes spontaneously, until the vibrating bases wash out the cluster. 

Moreover, we also observed some spontaneous tendency of the system to segregate the species, a fact which becomes 
more apparent at high densities. Both phenomena have their origin in the inelasticity. A possible qualitative expla- 
nation of the observed dynamics is as follows: particles with smaller restitution coefficient tend to group together, 
since the more energy they dissipate through 2-2 collisions, the denser the segregated domain becomes. Particles with 
higher restitution coefficient bounce for a longer time after a collision, and dilate more quickly. 



D. Velocity Probability Distributions 

Figure ^ shows the distributions functions for the translational velocities along the horizontal direction. For sake of 
comparison, all the velocities were rescaled by their mean square values. We notice that the two transversal velocity 
distribution functions deviate from a Gaussian and are fitted by /(c) = A/ (exp((3c a ) + exp(—{3c a )), with: A = 1.068, 
[3 = 1.74, a = 0.918 . 

The exponential tails appear to have a smaller slope than the theoretically predicted 3/2 value, for a uniform system 
stochastically driven ^0|. We notice the phenomenon already reported in |f[2[ : the rescaled velocity distribution of 
lighter species display higher tails. 

Figure [t] also shows the angular velocity distributions (rescaled by their mean square value, see caption). 

Our simulations show that such distributions have pronounced non-Gaussian tails. The angular velocity distribu- 
tions can be described by the same scaling function we used for translational velocity distributions. 
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FIG. 8: Probability distribution functions of the collision times for the two components (system with rn = 0.93, rn = 0.7, 
ri2 = 0.8). For comparison we plotted the corresponding distributions for an elastic system. The inset shows the same curves 
in bilogarithmic scale. 



E. Collision time distribution 



Recently Blair and Kudrolli utilizing high speed digital photography measured the collision statistics of grains 
bound to move on an inclined plane |2l| . They determined the distribution of path lengths, P(l) and showed that it 
deviates from the theoretical prediction for elastic hard spheres. In particular, P(l) shows a peak in the small I region, 
not present in elastic systems. In order to assess the existence of such a behavior in the IHS model we performed 
similar measurements. Figure || displays the results for the collision times of the two species. One clearly sees that 
the probability density that a particle suffers a collision in a short interval is enhanced with respect to the elastic 
case. The physical reason for that lies in the existence of strong correlations which lead to the presence of clusters 
where the path lengths are shorter than in a uniform system. For the sake of comparison we plotted in figure ^ also 
the corresponding distribution of an elastic gas having the same density and granular temperature of our inelastic 
mixture. The effect of the shorter average collision time is clearly visible. We also notice that the ratio of the collision 
frequencies is approximately equal to the the ratio of the average velocities, as expected from elementary kinetic 
arguments. 



IV. LONG SYSTEM 



In the case of a box of dimensions L x = 48 d and L y = 64 d with 300 grains of each species the effect of gravity 
is much more evident, resulting in a stronger inhomogeneity and asymmetry of the system. In fact, the density 
and temperatures profiles are not symmetric with respect to the vertical direction. A configuration of such system 
is shown in fig. ^. Most of the particles remain suspended above the bottom wall and are hit by those which are 
between the bottom wall and and the bulk. Very few particles reach the upper vibrating wall, so that the granular 
temperature of the system is much lower in the top than in the bottom. The partial temperature profiles and density 
profiles are shown in fig. |l0|. One sees that near the lower vibrating wall the temperature profiles are similar to those 
of the shorter system, whereas in the bulk they are appreciably different. The average ratio is lower than in the 
L y = 32 d case. An interesting feature present in fig.J^ is the presence of a region of increasing temperatures. Such 
a phenomenon has been predicted by the hydrodynamic theory of Brey ct al. H. The theory predicts a temperature 
varying as T ~ z 3 / 2 in the region above the minimum. Such a prediction is verified by our simulation results. 

Physically, the increase is due to the competition between the energy input from the walls and the small dissipation 
due to the reduced number of collisions associated with the low density region. This can be appreciated in the inset 
of figure [II] where the quantity £(z) = n(z)T(z) 3 ^ 2 is shown: £(z) is proportional to the energy dissipation due to 
inelastic collisions among particles, being the collision rate cx nT 1 ^ 2 and the average dissipated energy in a single 
collision oc T. The graph of £(z) reveals that the energy dissipation is much more relevant in the bottom region than 
in the middle and upper regions. 
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FIG. 9: Snapshot of the system described in IV. Open circles indicate particles of species 1, black circles particles of species 2. 
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FIG. 10: Temperature and area fraction profiles for the system discussed in section IV. In the inset the dissipation rate f (z) is 
shown. In the figure we compared our result against the hydrodynamic prediction T ~ z 3//2 (solid lines). The vertical dotted 
lines correspond to the rest position of the horizontal boundaries. 



Correspondingly the velocity PDF shows an interesting behavior (see figure [IT]). In fact, the shape of the rescaled 
PDF becomes narrower with the height. In the central region the measured PDF resembles the exponential PDF 
measured in the short system (figure 0), whereas in the bottom region the PDF has more extended tails. Such a lack of 
universality in the velocity PDf was noted experimentally by Blair and Kudrolli (22) . It was also noticed in : there 
it was shown that the tails of the PDF became broader when the dissipation rate was increased. Here the mechanism 
is similar: the broader PDFs are those measured in the bottom region, where the dissipation rate is higher. 
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V. CONCLUSION 



Summarizing, we studied a systems of inelastic particles in vertically vibrated containers and subject to the grav- 
itational field. The system was numerically investigated by using an Event-Driven dynamics, which affords the 
exploration of a wide range of dynamical parameters. We examined two similar setups which differ only for the aspect 
ratio of the container and for the number of particles. 

In the shorter system we determined the partial granular temperatures, their ratio, and the area fraction profiles 
along the vertical direction. Those measures are in qualitative agreement with the experimental results, but a 
quantitative comparison requires a more detailed knowledge of the experimental parameters. Physically the lack of 
energy equipartition between the two species has two causes: a) the vibrating walls feed energy proportionally to the 
mass of each species; b) the dissipation is in general different for the two components. The two effects may conspire 
(for example if the heavier species is the more elastic) to give a small value of the temperature ratio or they may cancel 
and give a temperature ratio close to one (for example when the heavier species is the more inelastic). Moreover, 
we observed that density profiles are non-uniform along the horizontal direction, as well, indicating that the particles 
tend to clusterize in the vicinity of side walls. We also measured the velocity distributions, verifying that they can 
be collapsed onto each other, under proper rescaling, and the scaling function has a stretched exponential behavior. 
Finally the distribution of flight times between succesive collisions has been measured and compared to that of an 
elastic system: the inelasticity has the effect of enhancing the statistics of very short times. 

In the longer system we have again obtained the temperature and area fraction profiles, observing a stronger inho- 
mogeneity and asymmetry. In this case the area fraction is much larger near the bottom wall, reaching higher values 
than in the previous experiment. The middle-upper region seems to be qualitatively well described by recent hydro- 
dynamic theories developed for one component systems. Here the velocity PDFs display a non-universal behavior, 
with broader tails in the more dense (and dissipative) regions. 



VI. APPENDIX 



In order to make it simpler for the reader to interpret the present model, we present an appendix in which we 
explicitly state the collision rules p^]. 

The colliding particles are characterized by radii R\ and i? 2 , positions i"i and r 2 , translational velocities vi and 
V2 and rotational (angular) velocities u}\ and w 2 (we assume that if oj is parallel and in the direction of the z axis, 
than the rotation is anticlockwise if seen from above the xy plane). We introduce the normal unitary vector joining 
the centers of the particles n = (r 2 — ri)/|r 2 — ri| and the tangential unitary vector t obtained rotating n by an 
anti-clockwise angle tt/2. Then we introduce the relative velocity g, the velocity of the center of mass V and the 
velocities of the particles in the center of mass frame Ci arL d C 2 : 
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g = vi - v 2 (2a) 

V = miVl+m2Vz (2b) 

mi + m2 v ; 

C. = (2c) 



mi 
m 2 



C 2 = -^g (2d) 



where m e // = mim^/ (rn\ +m 2 ). 

Then we decompose the relative velocity g on the orthonormal basis given n and t, as well as the velocities of the 
particles in the center of mass frame, i.e.: 



9n = (vi - v 2 ) • n (3a) 
9t = (vi - v 2 ) • t (3b) 

Cnf, = ^9n (3C) 



m eff 



Ctp = ^lLg t (3d) 



with j3 = 1, 2 the index of the particle. 

We finally introduce g c as the relative circular velocity at the point of contact and g r as the total tangential relative 
velocity (circular and translational) at the point of contact: 

g c = Riuji + i? 2 w 2 (4a) 
9r = 9c + 9t (4b) 

To characterize the collision rules we use a model that take into account a reduction of normal relative velocity 
(<?„), a reduction of total tangential relative velocity (g r ) and an exchange of energy between those two degrees of 
freedom. The reduction of normal relative velocity is modeled as usual by means of a restitution coefficient a € [0, 1]: 

9n = -ag n - (5) 
We assume a dependence of a by the relative velocity of the form: 



a{.9n) 



per g n > v 



(i-ao) 1 ^ 4 per g n < v 



where ccq is a constant which stands for r%j , vq — y gD, with D the average diameter of the particles (and g is the 
gravity acceleration). From equations (Q) we obtain the update of normal velocities in the center of mass frame: 

Ci = — —9n (6a) 
mi 

Cn2 = " —9n 

m 2 

What lacks now is an expression for the tangential and angular velocities after the collisions. We distinguish between 
two possible cases: sliding or sticking collisions. The condition that allows to determine if a collision is sticking or 
sliding is the following: 



(7) 
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where I is the adimensionalized inertia moment and is equal to 1/2 or 2/5 if the particle is a disk or a sphere 
respectively, while fi is the static friction coefficient of the surface of the particles which in the following will be 
assumed to be equal to the dynamic friction coefficient. 

In the sliding case we use the following rules to update the tangential components of the velocities of the particles 
in the center of mass frame: 



Cti = Cti - M(l + a)—^-g n sign(g r ) 
mi 

7Th f f 

C*2 = Ct2 + Ml + a )—^—9nsign(g r ) 

m 2 

, + a) m e tf 
R\uj 1 = Riloi : '-^-g-nsignygr) 

c / D At(l + a) m eff . 
R 2 uj 2 = R 2^2 : g n sign(g r ) 

I m 2 

In the case of a sticking collision, instead, the update rules are obtained considering that: 

Cti ~ Ct2 + Riu'i + R^2 = o 

from which, after calculations, one gets: 





1 

= Z + l' 


Ct2 


1 


= Z + l' 




= Riwl 




— R 2 ui 2 - 


R 2 u> 2 


= R 2 u2 




- Riui- 



Oi - 



m eff 



I + 1 mi 



(RlLUl + R 2 u> 2 ) 



I Tip f f , „ „ 

Ct2 + r— (Riwi + R2U2) 

l + 1 m 2 

I meff 



l + l (Z + l)m 2 

m eff 1 r 

- Ol 



(Z + l)mi l + l 
I m eff 



l + l {l + l)mi 



l eff 



1 



(/ + l)m 2 



+ 1 
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(8a) 
(8b) 

(8c) 

(8d) 

(9) 

(10a) 
(10b) 
(10c) 
(lOd) 
(lOe) 
(lOf) 



The velocity of particles in the absolute frame are finally obtained in the two cases (considering that the center of 
mass is not perturbed by the collision) by the equation: 



(for particle of index /?) leading to the following global collision rule for translational velocities: 



(11) 



vi = vi - (1 + a) 



v 2 = v 2 + (1 + a) 



m 2 



mi + m 2 



mi 



mi + m 2 



[(vi - v 2 ) • n]n • 



[(vi - v 2 ) • n]n- 



8ign(g r Ml + a )^T 2 \^ ~ ^2) ■ n]t (sliding) 

' ; ! ; - ) [(v 1 - V2 )-t]t- T ^ T ^-(R 1 u Jl +R 2 oj 2 ) (stick) 



(i+l)(mi+m 2 ) I 



(12a) 



sign{g r )n{\ + a) m ™ m2 [(vi - v 2 ) • n]t (sliding) 
jF^jfeKvi - v 2 ) • t]t + ^^(R^ + R 2 u 2 ) (stick) 

(12b) 



while for rotational velocities: 
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- ^ m r+m 2 9nsign(g r ) (stick) 

Z+l (l+l)(mi+m 2 ) 



R2W2 



R 2 uj 2 - M(1 + a) mi 



mi+m 2 
m2 



(i+l)(mi+m 2 ) 



(i+l)(mi+m 2 ) 

g n sign(g r ) (stick) 

p , . 2211 





^■Cti (slide) 
TT^T n,,,.. ,,,, 7 t j-prCn (slide) 



(13a) 
(13b) 
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